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Abstract A Bayesian method for forecasting solar cycles is presented. The 
approach combines a Fokker-Planck description of short-timescale (daily) fluc- 
tuations in sunspot number (No ble and W heatland, 2011, Astrophys. J. 732, 5) 
with information from other sources, such as precursor and/or dynamo models. 
The forecasting is illustrated in application to two historical cycles (cycles 19 
and 20), and then to the current solar cycle (cycle 24). The new method allows 
the prediction of quantiles, i.e. the probability that the sunspot number falls 
outside large or small bounds at a given future time. It also permits Monte 
Carlo simulations to identify the expected size and timing of the peak daily 
sunspot number, as well as the smoothed sunspot number for a cycle. These 
simulations show how the large variance in daily sunspot number determines 
the actual reliability of any forecast of the smoothed maximum of a cycle. For 
cycle 24 we forecast a maximum daily sunspot number of 166 ± 24, to occur in 
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March 2013, and a maximum value of the smoothed sunspot number of 66 ± 5, 
indicating a very small solar cycle. 

Keywords: Solar Cycle, Models; Sunspots, Statistics 
1. Introduction 

The solar cycle is the semi-periodic change in the level of magnetic activity 
on the Sun, driven by a cyclical change in the Sun's local magnetic field. The 
International Sunspot Numbciy, defined as 

s = k(10g + N), (1) 

where g is the number of sunspot groups, N is the number of individual spots, 
and k is a correction factor, is the most commonly used measure of solar ac- 
tivity ( |Bruzek and Durrant, 1977[ ). Variation in the sunspot number is com- 
prised of the semi-regular 11 year cycle, as well as large daily, weekly, and 
yearly fluctuations on top of the underlying secular variation QAbreu et al, 2 008; 
Noble and Wheatland, 20TT |. 

The sunspot number is an important indicator of solar activity, and hence 
of the space weather experienced on the Earth ( |Petrovay, 2010[ ). As a result 
reliable prediction of the sunspot number is important. There are several meth- 
ods for forecasting sunspot numbers (see |Petrovay| 120101 IKanel 120071 IPesnelll 
2008] for reviews), including precursor methods, time-series methods, and dy- 
namo model based methods. Precursor methods correlate aspects of a future 
solar cycle {e.g. sunspot maximum) with indices of current solar and geomag- 
netic activity. Time-series methods extrapolate sunspot data into the future 
using mathematical or statistical techniques, including nonlinear models {e.g. 

2 Sunspot data are provided by the US National Geophysical Data Center (NGDC) at 
|http://www . ngde. noaa.gov /stp /spacewea ther.html | 
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and neural networks (e.g. |Conway[ I1998|) . Physical predictions are provided by 
dynamo-based models, which start with models of the Sun's internal dynamo, 
the source of the magnetic fields of sunspots (e.g. Dikpati and Gilman| [2006; 



|Dikpati, de Toma, and Gilman| I5006P . Precursor, time series, and dynamo-based| 

methods all forecast solar activity with some success ([Hathawa y, 2009 1, with 



precursor methods being the most successful QKane, 2007[ ). 

Two specific criticisms of existing approaches to prediction are that it is 
difficult to rigorously combine forecasts from the competing methods, and that it 
is unclear how to update and/or reconcile forecasts (in particular long-range pre- 
cursor forecasts) with new sunspot data as it becomes available. These problems 
were addressed by |Hathaway, Wilson, and Reichmann] ()1999[) , who used precur- 
sor methods to forecast the size of the underlying cycle and then regressions 
to update precursor forecasts as new data became available. In this paper we 
also present an approach to prediction which combines precursor and time-series 
methods. Our method is similar to that of [Hathawa y, Wilson, and Reichmann 



(1999]), with some specific differences. We combine sunspot data with precursor 
forecasts in a statistically rigorous way using a Bayesian framework, rather than 
combining methods using weighted averages. We also forecast the short time- 
scale fluctuations in sunspot number (e.g. the variance in the daily sunspot 
number), in addition to the size and shape of the underlying solar cycle. 

Section [5J covers the theory of the Bayesian approach to forecasting daily 
sunspot numbers. Section 12.11 introduces the Fokker-Planck model used to de- 
scribe the distribution of sunspot numbers, and Section 12.21 gives the details of 
the Bayesian method. Section 12.31 discusses the historical average solar cycle, 
or 'mean cycle', which we use as a starting point for prediction. Section 12.41 
introduces a Monte-Carlo approach to simulating daily sunspot numbers based 
on analytic approximations to the Fokker-Planck model, which permits analysis 
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of the size and shape of the average underlying solar cycle, as well as construc- 
tion of the joint distribution of the size and timing of maximum daily sunspot 
number. We also show how the large variance in daily sunspot number results in 
large variance in the monthly smoothed maximum sunspot number (i?)max, and 
discuss the implications of this for the reliability of any forecast of the maximum 
of the cycle. Section [3] illustrates the Baycsian framework from Section [5] in 
application to two historical solar cycles, namely cycle 19 (in Section 13. 1|) and 
cycle 20 (in Section |3~2")) . Section 2] applies the Bayesian method to forecast the 
current cycle (cycle 24). Techniques from Section [2] are used to quantify the size 
of large/small sunspot numbers during cycle 24, and to estimate the likely size 
and timing of the next maximum in both daily sunspot number, and monthly 
smoothed sunspot number (i?)max- 

2. A Bayesian Approach to Solar Cycle Forecasting 
2.1. Fokker-Planck Model for Sunspot Number 

The solar cycle variation in sunspot number comprises long-term secular vari- 
ation, and large short-term statistical fluctuations. The long-term variation 
may be considered the underlying solar cycle, driven by the internal dynamo, 
and the short-term fluctuations attributed to complicated physical processes on 
the solar surface associated with sunspot formation, evolution and dispersion 
( |Parker, 19551 ). 

The long-term solar cycle variation over a single cycle may be described by a 
cycle amplitude, cycle period, and cycle asymmetry (Hathaway, Wils on, and Reichmann, 1994[ ),[ 
which we represent with a 'driver function' 8(t). For illustrative purposes we 
consider first the simple choice for 6(t) of harmonic variation: 

6{t) =a Q + a 1 sm{2Trt/a 2 + a 3 ), (2) 
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where a\ is the cycle amplitude, a-j « 11 years is the cycle period and 03 is the 
cycle phase. With this choice there is no cycle asymmetry. 

Short-term fluctuations in sunspot number on top of the driver function 6(t) 
may be modelled using a probability distribution function f(s, t) = f(s, t\so, to), 
such that f(s,t)ds is the probability that the sunspot number is between s and 
s + ds at time t. This approach represents the sunspot number as a continuous 
random variable. [Noble and W heatland! (|201ip modelled the time evolution of 
f{s,t) using the Fokker-Planck equation: 

8f 1 B 2 B 

~bl = 2d? [^*)fM] - IMM)/(M)] , (3) 

where fi(s,i) describes deterministic changes in sunspot number, and a 2 (s,t) is 
a variance describing stochastic variation. Sunspot number is non-negative, so 
the appropriate boundary condition at s = is a 'zero probability flux' condition 



H(s,t)f(s,t) - \§- s [<r 2 (s,t)f(s,t)] 



= 0. (4) 

s=0 



An appropriate choice for //(s,i), in terms of the driver function 9(t) is 

MM) =*[*(*)-*]■ (5) 

This choice causes the fluctuating sunspot numbers to tend to return to the 
value 6(t) with a characteristic timescale 1/k. In this way the driver function 
9(t) represents the secular or long-term sunspot number. 

The size of the observed squared deviations r 2 {t) = [s(t) — s(t — At)] 2 , a 
proxy for daily sunspot number variance, tends to increase with sunspot number 
( jNoble and Wheatland, 2011| ). Therefore a simple choice for the variance, which 
models this increase is 

cr 2 (M) = A) + /3is + /3 2 s 2 , (6) 
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where /3o describes variance in sunspot number at s = 0, and /3i and /?2 describe 
the increase in variance with sunspot number (we assume that /?o, Pi and /?2 are 
all positive or zero). The model then has four parameters (k, f3o, ft\ P2), together 
with any parameters in the driver function 9(t). 

The choices of Equations (j4]), ([5]) and (|6|) arc discussed in detail in lNoble and Wheatland j 
(20lT|). The authors showed that the model given by Equations @ to © applied 
to historical monthly sunspot number data generates a probability distribution 
function f(s, t) which agrees both quantitatively and qualitatively with observed 
sunspot statistics, even for the simple harmonic choice for 8(t). The model pa- 
rameters were estimated from the historical data using a maximum likelihood 
procedure explained in Section 12.21 

A more realistic choice of driver function 9(t) for a single solar cycle than 
Equation @ is provided by the functional form QHathaway, Wilson, and Reichmann, 1994[ ):| 



r ° (t ~ f0 2 )3 , 1 



exp 



where to is the start of the cycle, and a, b and c represent the cycle amplitude, 
period, and asymmetry respectively. With this choice of driver function there 
are seven parameters in the Fokkcr-Planck model, which may be represented in 
a vector 

n = [a,b,c,K,Po,Pi,fa]. (8) 

The distribution of model sunspot numbers is written f(s,t;fl) to indicate the 
explicit dependence of the distribution on the model parameters. We assume that 
the cycle start date to is known, but it could be treated as another parameter 
and estimated from the data. 

If the parameters ft generating sunspot data arc known, the time evolution 
of the distribution of sunspot numbers / (s, t; ft) is uniquely determined by the 
Fokker-Planck Equation (|3]). However, the parameters are unknown. To describe 
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historical sunspot numbers the parameters may be estimated from historical 
sunspot data, following INoble and W heatland (2011]) [who used Equation @ as 
the choice of 9(t)]. For forecasting, it is necessary to estimate values ft of the 
model to forecast future sunspot numbers. These procedures are explained in 
Section l2~2l 

2.2. Bayesian Estimation of Model Parameters 

Given an observed set s = {so, si, . . . , St} of sunspot numbers at times {to, ti, . . . , tr}M 
the maximum likelihood (ML) estimate is the parameter set ft which maximises 
the likelihood function representing the probability of the data given the model: 

i=T 

£(s\fl) = Ylf(si,U\si-uU-i;n). (9) 

We assume that the distribution of Sj at time ti depends only on the previous ob- 
servation Si-i at time tj_i, which is the Markov property ( |Karatzas and Shreve, 1988] ) .| 
ML estimates are optimal in the sense that they are both efficient and consistent 
in large samples ( |Dacunha-Castelle and Florens-Zmirou, 1986[ ). However, ML 
estimates are limited in that they only use information from the observed data, 
and ignore other information which may be available. 

With sunspot data we have additional information about the possible size, 
shape, and length of a future solar cycle, which may be included in the forecast 
using the Bayesian method (e.g. ISivial I2006p . The additional information may 
be in the form of a dynamo-based forecast, or a precursor forecast, for a future 
cycle. Our confidence in the reliability of this information is represented by a 
'prior distribution' V (ft) for the model parameters ft given the information. For 
example, if a precursor forecast for the variance parameter j3o is Poj an d if the 
parameter /3q is not correlated with other model parameters, then an appropriate 
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choice of a prior distribution for this parameter is 



V Wo) = 



2 ™l 



exp 



1 fPo-po 



2 V VPo 



(10) 



where the variance a\ represents how confident we are that /3q coincides with 



0o 



Because we are dealing with multiple parameters which may be correlated, 
it is necessary to include possible correlations in the prior. For the choice of 
Gaussian-distributed priors it is appropriate to use the general multinormal 
distribution 



1 



(27r) fc/2 |E|V2 



exp 



(n-n) s _1 (n-n)' 



(ii) 



where k is the number of parameters in f2, and the matrix E is the variance- 
covariance matrix describing the uncertainties in each parameter and the co- 
dependence of the parameters. For the parameters in thclHathaway, Wils on7and Reichmannj 
(|1994p driver function (Equation ([7])), E is a 7 x 7 matrix of the form 



'b,a 



\ 



(12) 



where <rf is the uncertainty in parameter i and <Tij is the covariance between 
parameters i and j, for i, j = a, b, . . . , 02- As an example of the importance of 
correlations, it is well known that cycles which rise rapidly tend to be large 



(Wa ldmeier, 1935 1, so that we expect parameters describing the period and 
amplitude to be negatively correlated. 

Predictions incorporating prior information may then be made as follows. 
When a cycle begins, daily sunspot data s = {sq, s±, . . . , st} becomes available. 
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This can be combined with the prior information by calculating the 'posterior 
distribution' V (O | s) , according to Bayes' rule QSivia, 2006[ ): 



p(«|s) 



r(s\n)v(n) 



(13) 



The term V (s|fi) in Equation (fTU|) is the likelihood function ©, and the denomi- 
nator is a normalising constant. The posterior distribution combines information 
about l~i contained in the data (the time series approach), with relevant infor- 
mation external to the data (e.g. from precursor and/or dynamo models, or any 
other source). 

We are unable to solve the Fokker-Planck equation © analytically, so we 
cannot evaluate the likelihood V (s[J2) in closed form. An analytic approximation 
appropriate for daily data is ( |Noble and Wheatland, 2011 ): 



f(s,T\s ;ft) 



1 



^2ir<T 2 (s ,t )i 



exp • 



(s + /j,(s ,to)r) 
2a 2 (s ,t )T 



- exp ■ 



[S+(S Q + fJ,(s ,t )T)Y 

2a 2 (s ,t )T 



(14) 



where r = t — to- Equation (|14[) is the conditional probability distribution func- 
tion of the random variable where s(t) is a normal random variable with 
mean sq + /i(so,io) and variance a 2 (so,to)T, and sq is the sunspot number at 
time to. Using this approximate solution, the likelihood function © is 



± 

V{s\Sl) = (2ny^l 



exp 



[Si - (Si_i + Mi-l r ) 



2a 



■ exp ■ 



[Si + (Sj_i + (H-lT 



i-l • 

2 1 



2a 



i-i' 



(15) 



where /ij = /i(sj,£j) and cr? = a 2 (si,ti). 

Given the posterior distribution V (fl\s) a specific estimate for the parameters 
fZ may be calculated in a number of ways flJaynes and Bretthorst, 2003] ). In this 
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paper we use the the most probable, or 'modal estimate' of SI: 

ft = argmaxT? (Sl\s) . (16) 
2.3. Construction of a Mean Solar Cycle Prior 

Rather than using a specific precursor forecast, we consider using an historical 
average solar cycle, which we refer to as a 'mean cycle', as a prior. This means 
that before the start of a cycle the most probable shape of the cycle (i.e. the 
parameters in the driver function) and the variance of the cycle (i.e. the param- 
eters in a 2 (s,t) of the coming cycle) are represented by an historical average. 
This choice may be interpreted as a 'guess in total ignorance'. 

To determine the parameters f2 for the mean cycle we consider daily sunspot 
data for the previous 13 solar cycles over the interval 1850 to 2010. The Hathaway, Wils on, and Reichmann| 
(1994]) driver function given by Equation ([7]) is assumed to represent the shape 
of each underlying cycle, and the variance of each cycle is modelled by Equation 
([5]). For each solar cycle maximum likelihood estimates f2 = a, b, c, k, /?o, fti, $2 
of the seven model parameters are calculated, as shown in Table [TJ 

The average for each parameter over the previous 13 cycles is assumed to 
represent the mean cycle, and is used in our prior distribution. The sample 
means (denoted Cl) are given in Tabled The variance-covariance matrix E is 
calculated using sample covariances between the ML estimates for the seven 
parameters in Table [1] For example the covariance between the amplitude a and 
period b in £ is 

13 

12 



< J a,b = -^Yl(ai-a)(bi-b). (17) 

i=i 

The (non-dimensional) correlation matrix is given in Table |31 A number of 
important correlations exist. In particular, the large correlation between k and 
Pi (89%) shows the strong relationship between the variance parameters and the 
rate at which sunspot number returns to the level 9(t). There are also significant 



SOLA: SP_3.tex; 15 November 2011; 1:19; p. 10 



A Bayesian Approach to Forecasting Solar Cycles Using a Fokker-Planck Equation 



correlations between the size of the cycle a, and the time to maximum b and 



asymmetry c (collectively describing the Waldmcicr Effect ( Waldmeier, 1935 1) 
2.4. Fokker-Planck Modelling of the Mean Solar Cycle 



In this section we investigate characteristics of the average solar cycle using the 
mean cycle parameters Jl estimated from daily sunspot data over the interval 
1850 to 2010, given in Table [2] The driver function Equation (JT]) with parameter 
values from Table [2] describes the average size and shape of the underlying solar 
cycle. Short-term deviations in sunspot number from this average are described 
by the three variance parameters in Table O 

Figure Q] illustrates the mean cycle (i.e. the sunspot model with f2 = fi), 
and simulations of daily sunspot number over the mean cycle using the Fokker- 
Planck model. The red curve (solid) in Figure[T]is the |Hathaway, Wilson, and Reic hmann|| 
(199||) driver function 0(t), given by Equation (0 with mean cycle parameters 
a, b and c. The maximum value of the driver function $max = 119 occurs 4.4 
years after the cycle start date. To investigate the likely size of short-term devia- 
tions about the driver 0(t), we numerically solve the Fokker-Planck equation ([3]) 
for the initial initial condition f(so,to] = S(so) with s(io) = at to = 0. The 
solution is obtained from t = to t = 11 years. Based on the numerical solutions 
for f(s, t\so, to] Cl) we calculate the upper and lower 1% quantiles, which are the 
curves su(t) and Si(t) defined by 

Si(t) roc 

ds'f (s',t\s , t ; n) = / ds'f (s',t\s , to; n) = 0.01. (18) 

J su{t) 

These curves delineate boundaries of extreme sunspot numbers (i.e. the probabil- 
ity that the sunspot number is larger or smaller than sjj(t) or s^(t) respectively, 
from a given initial condition is 1%). The blue curves (dot/dashed) in Figure 
[1] show these upper and lower 1% quantiles for the mean cycle. The maximum 
value attained by the upper 1% quantile is su(t) = 234, which occurs 4.3 years 
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after the start of the cycle. This means that on average solar maximum occurs 

4.3 years after the start of the cycle, and that there is a 1% chance of observing 
a daily sunspot number larger than 234 at solar maximum for an average cycle, 
given s(i ) = 0. 

Using the analytic approximation (Equation ([14])) wc can simulate daily 
sunspot numbers over the mean cycle. To do this we repeatedly generate random 
variables s(t) from the conditional probability distribution (Equation (j 14|) with 
parameter values from Table [5J These numbers represent a sequence of possi- 
ble daily sunspot numbers (one Monte Carlo simulation). The green points in 
Figure [1] are an example of simulated daily sunspot numbers. The maximum 
daily sunspot number s* for this particular simulation is s* = 266, which 
occurs 4.4 years after the cycle begins. In this simulation 1.1% and 1.8% of 
the sunspot numbers are greater than and less than the upper and lower 1% 
quantiles respectively. 

With the Fokkcr-Planck model we can investigate the likely size and timing of 
daily sunspot maximum using repeated Monte-Carlo simulations. Wc denote by 
t* the time of the occurrence of the maximum s* of the daily sunspot number, for 
one simulation. Figure [2] shows a Monte Carlo estimate of the joint distribution 
f{s*, t*\so, to; Si) of the size and timing of daily sunspot maximum based on 
5 x 10 5 simulations. Each simulation is generated in the same way as the single 
instance shown in Figure [TJ The expected size of daily sunspot maximum (the 
average over the simulations) is (s*) = 271, which occurs approximately (t*) = 

4.4 years after the start of the cycle. This is comparable to the sample average 
maximum sunspot number from the previous 13 cycles (see Section [2. 3[) . which 
is 255. The largest value of the daily sunspot number of the 5 x 10 5 simulations is 
s* = 504, which suggests that mean cycle can generate extremely large sunspot 
numbers, although it is very unlikely (the mean cycle is expected to generate 
one such event every 5 x 10 5 cycles, or 5 million years). 
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We can also investigate the monthly smoothed sunspot number (i?)max, 
which is the main focus of much of the existing literature ( |Petrovay, 2010[ ). 
During each simulation of daily sunspot numbers discussed above, we calcu- 
late a 13-month boxcar average sunspot number (R), and store the size of the 
maximum (-R)max- Figure [3] plots the distribution of (-R)max based on 5 x 10 5 
simulations of daily sunspot number. The expected value over the simulations is 
(-R)max = 125 ± 8, and the lower and upper 5% quantiles are are 113 and 138 
respectively. For comparison the average smoothed maximum from the previous 
13 cycles is 121. 

Figure [3] has important implications for any forecast of (-R)max- Even if the 
model parameters for a solar cycle are known, large daily variation in sunspot 
number causes large variation in the possible smoothed maximum value of the 
cycle. This indicates that the reliability of any forecast of (-R)max is limited by 
the large daily fluctuations in sunspot number. 

3. Forecasting Historical Solar Cycles 

In the following sections we consider application of the model to forecasting two 
historical solar cycles: cycles 19 and 20. Cycle 19 is chosen because it is the 
largest solar cycle since daily sunspot number records began QKane, 2002[ ). In 
contrast, Cycle 20 is very similar in amplitude and shape to the mean cycle. As 
such these cycles are very different in size and shape, and useful as illustrations 
of the Bayesian forecasting method described in Section 12.21 The two cycles are 
shown in Figure HI The maximum of the observed smoothed sunspot number for 
cycle 19 is 201, and the maximum for cycle 20 is 110. 

In the forecasts in this section we use the mean cycle constructed in Section 
12.31 as a prior, and then we make updated predictions using successively more 
historical data from the start of a cycle, to demonstrate the Bayesian prediction 
method from Section l2~2l We use the analytic approximation (Equation to 
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the solution to the Fokker-Planck equation to evaluate the likelihood function 
(|15[) . and we use the mean cycle estimates discussed in Section 12.31 with the 
multivariate normal prior distribution (Equation (fTTj) ). Because V(s) in Bayes' 
rule fj 13[) is required only as a normalising constant, we calculate the posterior 
distribution P(fi|s) oc 7- > (s|Jl)7 : '(f2), and evaluate the modal estimate fl of 
Equation (|16[) by numerical determination of the location of the maximum of 
the posterior distribution. 

3.1. Solar Cycle 19 

Solar cycle 19, which occurred from 1954 to 1965, is the largest cycle on record. 
As such it provides a useful illustration of how forecasts starting from a prior 
consisting of the mean cycle are modified by observation of larger than expected 
sunspot numbers. 

First we consider applying the Bayesian estimation procedure to sunspot data 
for the entire cycle. If we take the mean cycle as the prior and take all daily 
sunspot numbers from 1 January 1954 to 31 December 1964 as data s, construc- 
tion of the posterior V(£l\s) oc V(s\tt)P(n) and estimation of parameters gives 
the modal estimates ft in Table |4] The difference between the ML estimates 
in Table Q] and the Bayesian estimates in Table @] is the influence of the prior 
distribution in the calculation of the posterior. 

We can repeat the calculation using only part of the data from the start of 
cycle 19 as input in s. By constructing the posterior repeatedly, using successively 
more data from the start of the cycle, we mimic the process of forecasting and 
updating the forecast. Figure [5] illustrates the process of successive forecasting 
for this cycle. The green points are the observed daily sunspot numbers for 
cycle 19. The driver function for the mean cycle (the prior for the forecasts) is 
shown in blue (dot/dash). The three black curves (solid) are the driver functions 
given by Equation ([7]) calculated using Bayesian model parameters estimated 
with different amounts of daily sunspot data. The black curve with the smallest 
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maximum value is obtained using ten days of data from the start of the cycle, 
the next smallest uses one year of data from the start of the cycle, and the black 
curve with the largest maximum value uses two years of data. The driver function 
corresponding to the final Bayesian estimate using all data for the cycle (which 
has a maximum of #max = 182) is shown by the red dashed curve. This figure 
shows how, as parameter estimates are updated with additional daily sunspot 
data, the size, period and asymmetry of the forecast of the underlying solar 
cycle changes. Initial estimates of the size of the sunspot maximum are lower 
than that of the mean cycle, but the large sunspot numbers observed from about 
1956 onwards cause the estimates of the cycle maximum to increase. 

Figure |5] shows the estimate of maximum smoothed sunspot number (i?)max 
as a function of the time of the latest data used for the prediction. These 
estimates are calculated by averaging over 10 simulations (blue squares), as 
discussed in Section 12.41 The first estimate of (-R)max is calculated using 10 
days of data which consisted of 10 consecutive spotless days. The solid black 
curve is the expected value of (i?)max calculated using all daily data for cycle 
19. Early Bayesian estimates (i.e. using data from 1954 to 1955) of (-R)max are 
small because the data is dominated by a large number of days early in the cycle 
with zero sunspot number. From 1955 onwards the sunspot numbers increase 
more rapidly than expected for the mean cycle. As a result the Bayesian result 
for (i?)max rises rapidly until around 1958, and after that the estimate of (i?)max 
is approximately constant, fluctuating between 180 and 195. The final estimate 
(i.e. the estimate using all daily sunspot data for cycle 19) is (-R)max = 189 ±11, 
corresponding to the parameters in TableS] The observed value of (-R)max = 201 
is shown by the black dashed line, and is roughly one standard deviation higher 
than the expected value, given the data. This difference illustrates the large 
variability in the cycle maximum possible due to the daily sunspot number 
fluctuations (see Section |2~4")) . 
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3.2. Solar Cycle 20 

Solar cycle 20, which occurred from 1965 to 1976, is substantially different in 
character to cycle 19, discussed in Section 13.11 The shape of this cycle is more 
typical, similar to the mean cycle. 

Following the approach in Section 13.11 we first consider Bayesian estimation 
applied to daily sunspot data for the entire cycle, using the mean cycle as a 
prior. The data span 1 January 1965 to 31 December 1976. Table [5] lists the 
model estimates Q for the Fokker Planck model parameters obtained using the 
Bayesian procedure from Section 12.21 applied to the daily sunspot data for the 
whole cycle. The difference between the Bayesian and ML estimates is again due 
to the influence of the prior in the calculation of the posterior distribution. 

Figure [7] illustrates predictions for cycle 20 following Figure [SJ The green 
points are the observed daily sunspot numbers for cycle 20. The three black 
(solid) curves show the driver function (Equation (0) calculated using suc- 
cessively more sunspot data. The black curve with the largest maximum is 
estimated using ten days of data from the start of the cycle, the next largest uses 
one year of data, and the third black curve uses two years of data. The driver 
function corresponding to the final Bayesian estimate (which has a maximum 
#max = 124) is shown by the red dashed curve. Initial estimates of the cycle 
amplitude are larger than that of the mean cycle, as indicated in Figure [7] The 
observation of many days with small sunspot numbers (i.e. Sj < 50) up to three 
years into the cycle causes these large initial estimates of cycle amplitude to be 
reduced. The timing of the maximum of the cycle is correspondingly adjusted 
from late 1969 to late 1968. 

Figure [5] shows the estimate of maximum smoothed sunspot number (-R)max 
as a function of the time of the last data used for the prediction (blue squares) . 
These estimates are calculated by averaging over 10 5 simulations, as discussed 
in Section HOI The solid black curve is the expected value of (-R)max calculated 
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using all daily data for cycle 20. In this case there are a significant number of 
days with relatively large sunspot number at the start of the cycle (1965-1966), 
which cause the initial estimates of (-R)max to be larger than that of the mean 
cycle. However, from mid 1966 onwards there are many days with small sunspot 
numbers (i.e. s, < 50), and few days with large sunspot number (i.e. Sj > 0(£j). 
This causes the forecast to be reduced. The declining phase of cycle 20 (1969- 
1972) features a significant number of days with large sunspot number, which 
causes the forecast to increase again. The final estimate of (-R)max using all 
daily sunspot data for cycle 20 is (i?)max = 133 ± 11, corresponding to the 
parameters in Table[5]). The observed value of (i?)max is 113, which is roughly 
two standard deviations less than expected, given the data, again illustrating 
the possible variability in the maximum value. 



4. Forecasting the Current Solar Cycle (Cycle 24) 

In this section the Bayesian forecasting procedure is applied to forecasting the 
current solar cycle, cycle 24. There is considerable interest in forecasts for the new 
cycle given its late start and slow early onset ( Russell, Luhmann, and Jian, 2010] ) .[ 



In particular the years 2008 and 2009 featured long sequences of days in which the 



Sun was 'spotless' (Tokumaru et ai, 20091, and various features of the new cycle 



have prompted speculation that future activity will be substantially reduced (e.g. 
|Livingston and Penn] 12009ft 

Following Sections 12.21 12.31 and 12.41 the mean cycle is used as a prior. The 
start of cycle 24 is taken to be 1 January 2009. With these assumptions the 
Bayesian estimates of the Fokker-Planck model parameters using all available 
daily sunspot data 1 January to 31 March 2011 are given in Table [51 The 
maximum value #max of the driver function corresponding to the parameters 
in Table [6] is 61, which is approximately half the value for the mean cycle. The 
available data suggests that cycle 24 will be significantly smaller than average. 
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Figure |H] illustrates the forecasts for cycle 24 following Figures [5] and [7] The 
daily sunspot numbers for cycle 24 for the interval January 2009 to March 2011 
are shown by the green points, and the driver function for the mean cycle is shown 
by the blue dot-dashed curve. The solid black curve with the largest maximum 
value is the driver function using the Bayesian estimate of the Fokkcr-Planck 
model based on the first year of sunspot data (January 2009 - January 2010), 
and the second solid black curve uses the first two years of data. Combining the 
mean cycle prior with the first year of data gives estimates of the driver function 
very similar to the driver function of the mean cycle. However, due to the large 
number of days during the latter part of 2010 with small sunspot numbers, the 
driver function using the first two years of data has a much smaller maximum 
#max than that of the mean cycle. The red dashed curve is the forecast using 
all data, which has a maximum #max = 61. 

Figure [TU] shows the expected value of (i?)max as successively more data 
are incorporated into the Bayesian prediction method starting from 1 January 
2011, following Figures [5] and [5J These (blue squares) estimates are calculated by 
averaging over 10 5 simulations, as discussed in Section l2~4l The solid black curve 
is the expected maximum of (i?)max for the mean cycle. The early forecasts 
of (i?)max are lower than that of the mean cycle because of the significant 
number of days during 2009 with zero sunspot number. The forecasts of (i?)max 
steadily increase from early 2009 until mid-2010, but sunspot activity defied 
expectation and did not significantly increase during the latter part of 2010, and 
this causes a dramatic reduction in the forecast for (-R)rnax during late 2010 
and early 2011. The final estimate using all available data (and matching the 
parameters in Table |6]) is 2009 is (-R)max = 66 ± 5. This suggests that cycle 
24 will similar in size to cycle 14, and thus larger only than cycles 5 and 6. 
This prediction is close to the smaller estimates in the literature. For example, 
|Aguirre, Letellier, and Maquet| (|2008p predicted a smoothed sunspot maximum 
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of 65 ± 16, ICameron and S chtisslcrl (|2007p predicted a smoothed maximum of 
69 ± 15, and lKakadipTTTj) a smoothed maximum of 74 ± 10. 

Figure [TT] provides a third representation of the forecasts for cycle 24 based 
on the daily sunspot data for January 2009 to March 2011 (the observed data 
are shown in blue). The solid red curve is the driver function forecast based 
on all observed data, matching the Baycsian model estimates in Table |51 From 
April 2011 to January 2019 the solid red curve provides a basis for prediction of 
the upcoming sunspot numbers. The two dot-dashed black curves are the upper 
and lower 1% quantiles for the sunspot number distribution for the forecast, 
defined by Equation (JTSj) . These quantiles show the probability of excursions to 
large and small daily sunspot numbers. The maximum value attained by the 
upper 1% quantile is 138 during the period January-March 2013, which may be 
taken as the most likely time t* of daily sunspot maximum s*. The green points 
are simulated daily sunspot numbers for the remainder of cycle 24 using the 
Bayesian model estimates in Table |6l with initial condition s = 66 on March 31 
2011 (the sunspot number observed on that day). In this particular simulation 
the maximum daily sunspot number is s* = 168 which occurs during October 

2012. For the simulation a total of 0.9% and 1.1% of simulated sunspot numbers 
fall above and below the upper and lower 1% quantiles respectively. 

Figure [12] shows the joint distribution of the time t* and size s* of daily 
sunspot maximum for cycle 24, generated using 5 x 10 5 simulations of daily 
sunspot number based on the Bayesian estimates in Table [6] Averaging over the 
simulations we calculate the expected size of the maximum daily sunspot number 
to be (s*) = 166 ± 24, and this is expected to occur at a time (t*) during March 

2013. The sample average daily sunspot number maximum over the previous 13 
cycles is s* = 255, so on this basis cycle 24 is expected to be significantly smaller 
than average. The model probability that daily sunspot number maximum for 
cycle 24 is larger than the average maximum daily sunspot number s* = 255 is 
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V (s* > 255) = 0.4%. Hence it is unlikely that individual days with very large 
sunspot numbers will be observed during cycle 24. 

5. Discussion 

This paper introduces new techniques for estimating, analysing, and forecasting 
solar cycles, in particular daily and smoothed sunspot numbers for a cycle, 
and their statistical properties. In particular, we have shown that even with 
perfect knowledge of the details of a solar cycle, the observed sunspot maximum 
(either daily or smoothed) could achieve a broad range of values due to the large 
fluctuations in the daily sunspot number. This is important for all prediction 
done a priori, and indicates the true reliability of any forecast of the maximum 
of a cycle made before the fact. 

The main result of this paper is a new Bayesian prediction method for daily 
sunspot number fSection 12. 2|) . This method is illustrated in application to two 
dissimilar historical cycles (Section [3]) , and then is applied to the upcoming 
solar cycle (Section The method uses as a prior a mean cycle based on 
the observed solar cycles for 1850-2010 (Section 12. 3|) . Our investigation of this 
provides a characterisation of solar cycle variability which should also be useful 
to other workers. 

We model the sunspot number as a continuous-time stochastic process, with a 
probability distribution function described by a Fokker-Planck equation ( |Noble and Wheatland, 2011[ ) .| 
The Bayesian approach to forecasting uses the Fokker-Planck model to include 
information about solar cycles contained in sunspot data observed up to a 
given time, and combines these with external information (in principle that 
provided by precursor or dynamo-based forecasts). The external information 
is included by specifying an appropriate prior distribution. In this paper we 
take an historical average solar cycle ( a 'mean cycle') as a prior, which can be 
interpreted as a 'best guess in total ignorance'. However, the methodology can 
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accommodate any choice of prior. The Bayesian estimation method, combined 
with the Fokker-Planck equation approach, allows forecasts of the size and shape 
of the underlying solar cycle, as well as assigning probabilities to the observation 
of large deviations in sunspot number via calculation of upper and lower quantiles 
for future sunspot numbers. 

In addition, the Fokker-Planck model permits daily sunspot numbers to be 
simulated over a solar cycle, allowing Monte Carlo construction of the joint 
distribution of the size and timing of the maximum in daily sunspot number, 
as well as the distribution of the size of smoothed sunspot maximum (i?)max- 
In particular, the distribution of daily sunspot maximum determines the pos- 
sible size and timing of extreme sunspot numbers during a cycle, which define 
likely times for the occurrence of intense solar activity. Large flares and coronal 
mass ejections occuring at these times are drivers of our local space weather 
dCommittee On The Societal and Economic Impacts Of Severe Space Weather Events, 2008[ ),| 
and forecasting of extreme events space weather is an important task ( |Petrovay, 2010[ ).[ 

The application of the new method to the current solar cycle, cycle 24, 
provides insight into what we might expect over the next few years. Taking 
the mean solar cycle as prior and using data for 1 January 2009 to 31 March 
2011, the model forecast for the maximum of the smoothed sunspot number is 
(-R)max = 66 ±5, which is a very low value. The forecast maximum daily sunspot 
number is 166 ± 24, expected to occur during March 2013, and this is also very 
low. These predictions are consistent with other predictions in the literature in 
suggesting a much smaller than average cycle. The lack of a rapid rise in sunspot 
number during 2010, in particular, is shown by our modelling to imply a very 
small upcoming solar cycle. 
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Figure 1. Fokker-Planck modelling of the mean solar cycle, showing the driver function (solid 
red curve), upper and lower 1% quantiles blue dot-dashed curve), and an example simulation of 
sunspot numbers (green points), as described in Section 12.41 The maximum value attained by 
the upper 1% quantile is su(t) = 234, which occurs 4.3 years after the start of the cycle. In this 
simulation the maximum of the daily sunspot number is s* = 266, which occurs approximately 
4.4 years after the cycle begins. 

Waldmeier, M.: 1935, Astron. Mitt. Eidgen. Sternw. Zurich 14, 105. 
Yule, G.U.: 1927, Phil. Trans. Roy. Soc. A 226, 267. 
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Figure 2. The joint distribution of daily maximum sunspot number s* and time of maximum 
t* for the 5 X 10 5 Monte Carlo simulations of the mean cycle described in Section 12.41 The 
expected value of the maximum is 271, which occurs approximately 4.4 years after the start 
of the cycle. The largest daily maximum value in any simulation is 504, which suggests that 
the mean cycle has the potential to generate extremely large sunspot numbers, although it is 
very unlikely. 
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Figure 3. Distribution of the maximum of the 13 month smoothed sunspot number (-R)max 
calculated using 5 X 10 5 simulations of daily sunspot number for the 'mean cycle' (see Section 
12.311 , The expected value of the maximum is (i?)max = 125 ± 8. The upper and lower 5% 
quantilcs arc 113 and 138 respectively. 
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Figure 4. The daily sunspot number observations for cycles 19 and 20 (1954 to 1976) used 
for forecasting in Section [3] The red curve is a smoothed sunspot number. The maximum of 
the smoothed sunspot number for cycle 19 is 203, and the maximum for cycle 20 is 113. 
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Figure 5. Prediction of cycle 19 using successively more data from the start of the cycle. 
Daily sunspot numbers for the cycle are shown by the green points. The driver function for 
the mean cycle prior which uses no sunspot data is in blue (dot— dashed) , and the Bayesian 
estimate using all sunspot data for 1954 to 1964 is the red dashed curve. The model parameters 
for the red curve are given in Table [J] The solid black curve with the smallest maximum value 
is the forecast using ten days of data, the next smallest uses one year of data, and the largest 
solid black curve uses two years of data. 
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Figure 6. Prediction of cycle 19 using successively more data from the start of the cycle. The 
value of the maximum (-R)max = 203 for the observed cycle 19 data is shown by the dot-dashed 
line. The expected value of (-R)max calculated by average over 10 5 cycles is (-R)max = 189, 
and is shown by the solid line. The forecasts of (-R)max for the Bayesian modal estimate 
using daily sunspot up to the indicated time, and calculated by averaging over 10 3 simulations 
are shown by the squares. From 1955 to 1957 the forecast of (iJ)max rises rapidly and then is 
approximately constant. The final value, matching the parameters in TableQ] is (R)max = 182. 
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Figure 7. Prediction of cycle 20 using successively more data from the start of the cycle. Daily 
sunspot numbers for the cycle are shown in green. The driver function for the mean cycle prior 
which uses no sunspot data is the blue dot-dashed curve, and the Bayesian estimate using all 
sunspot data for 1965 to 1975 is shown by the red dashed curve. The model parameters for 
the red dashed curve are given in Table [5] The black solid curve with the largest amplitude is 
estimated using ten days of data, the next largest uses one year of data, and the largest solid 
black curve uses two years of data respectively. 
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Figure 8. Prediction of cycle 20 using successively more data from the start of the cycle. 
The value of the maximum (-R)max = 113 for the observed cycle 20 data is shown by the 
dot— dashed line. The expected value of (i?)max calculated by averaging over 10 5 simulations 
is (fl)max = 133, shown by the solid line. The forecasts of (i?)max for the Bayesian modal 
estimate using daily sunspot data up to the indicated time, and calculated by averaging over 
10 3 simulations arc shown by the squares. The significant number of large sunspot numbers 
at the start of the cycle cause early estimates of (iJ)max to be larger than expected. However, 
the lack of large sunspot numbers during solar maximum cause estimates of (i?)max to be 
reduced. The large variation in daily sunspot numbercauses estimates of (iJ)max to slowly rise 
to the final value ((ii)max = 133), matching the parameters in Table [5] 
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Figure 9. Prediction of cycle 24 using successively more data from the start of the cycle. 
Sunspot data for the cycle are shown by the green points. The driver for the mean cycle prior 
which uses no sunspot data is shown by the blue dot-dashed curve, and the driver function 
for the Bayesian estimates using all available data for the cycle (January 2009 to March 2011) 
is shown by the red dashed curve. The solid black curve with the largest maximum value is 
the forecast for the driver function using one year of daily sunspot data from the start of the 
cycle, and the second solid black curve uses two years of data. 



150 
140 
130 
120 

x 

J 110 

iSf 100 

90 
80 
70 



D09 



2009.5 



2010 

Year 



2010.5 



2011 



Figure 10. Prediction of cycle 24 using successively more data from the start of the cycle. The 
expected maximum (-R)max = 125 for the mean cycle is shown by the solid line. The forecasts 
of (fl)max for the Bayesian modal estimate using daily sunspot data up to the indicated time, 
and calculated by averaging over 10 3 simulations are shown by the squares. The lack of large 
sunspot numbers in late 2010 causes a dramatic reduction in the expected size of (-R)max- The 
final value, matching the parameter estimates in Table [6] is (i?)max = 66. 
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Figure 11. Prediction of cycle 24: illustration of the forecast distribution of daily sunspot 
numbers for the remainder of cycle 24. The model parameters used in the forecast are the 
Bayesian estimates given in Table [6] The solid red curve is the forecast of the driver function, 
and the two dot-dashed black curves are the upper and lower 1% quantiles for the sunspot 
number distribution. The blue points arc the daily sunspot numbers observed for January 2009 
to March 2011 used for the prediction. The green points are a simulation of future sunspot 
numbers using the parameters in Table [6] The upper quantile attains a maximum value of 138 
during the period January-March 2013, identifying this as the most likely time for a maximum 
in the daily sunspot numbers. 
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Figure 12. The joint distribution of sunspot maximum s* and time of maximum t* for a 
sequence of 5 X 10 Monte Carlo simulations of solar cycle 24, as described in Section [4] The 
simulation uses the Bayesian model parameters in Table [6] and the initial condition so = 66. 
The expected size of the daily sunspot maximum is 166 ± 24, which is most likely to occur 
around March 2013. 
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Table 1. Maximum likelihood estimates of the parameters CI = [a, b, c, k, /3rj, 01, 02] for thel 
previous 13 solar cycles over the interval 1850 to 2010. 



Cycle 


a 

[day 3 ] 




6 

[day] 


c 


K 

[day" 1 ] 


00 
[day" 1 ] 


01 

r i 1 i 

[day x ] 


02 
[day- 1 ] 




23 


7.8234x10" 


8 


1514.8 


0.22174 


0.08555 


17.689 


1.6569 


2.1862 x 10" 


5 


22 


14.112x10" 


8 


1368.2 


0.33153 


0.07305 


22.361 


1.6302 


1.1872 x 10" 


-3 


21 


12.497x10" 


8 


1414.0 


0.48998 


0.07325 


22.069 


1.4652 


3.0413 x 10" 


-3 


20 


7.2861x10" 


8 


1450.4 


0.90401 


0.06151 


44.288 


1.1235 


1.2022 x 10" 


-3 


19 


14.048x10" 


8 


1391.4 


0.66332 


0.07994 


18.636 


1.8863 


7.3281 x 10" 


-7 


18 


10.592x10" 


8 


1411.5 


0.65807 


0.09238 


22.541 


2.3351 


1.6617 x 10- 


-4 


IV 


7.2515x10" 


8 


1440.4 


0.80478 


0.11587 


30.236 


1.7644 


7.4863 x 10- 


.3 


16 


5.1614x10" 


8 


1457.8 


0.43823 


0.13072 


14.619 


2.8419 


4.0631 x 10" 


.3 


15 


7.8341x10" 


8 


1285.2 


0.83355 


0.10493 


23.245 


3.4815 


2.8562 x 10" 


.3 


li 


3.5130x10" 


8 


1576.8 


0.42673 


0.12617 


10.136 


3.5719 


8.5559 x 10" 


-4 


13 


8.5060x10" 


-8 


1256.3 


0.81241 


0.13579 


24.716 


3.3131 


2.4606 x 10" 


-4 


12 


3.7627x10" 


8 


1526.9 


0.55401 


0.12793 


14.543 


3.2136 


2.1382 x 10" 


-4 


11 


9.7135x10" 


8 


1356.6 


0.73723 


0.17509 


40.229 


4.5763 


9.1951 x 10" 


-4 



Table 2. Sample means fl for each model parameter estimated for the last 13 cycles. 
The solar cycle with SI = A is the mean solar cycle. 

a be ii 0o Pi 02 

[day- 3 ] [day] [day- 1 ] [day" 1 ] [day- 1 ] [day- 1 ] 

8.6233X10- 8 1419.3 0.60582 0.10633 23.486 2.5277 1.7123 xlO" 3 



Table 3. Correlation matrix for the model parameters estimated for the 
previous 13 solar cycles 1850-2010. 



a b c K fio 0i P2 



a 1.0000 -0.5268 -0.0243 -0.4671 0.2148 -0.4035 -0.1543 

b 1.0000 -0.5270 -0.0842 -0.3877 -0.1800 -0.0196 

c 1.0000 0.1596 0.6627 0.1721 0.2069 

re 1.0000 -0.0022 0.8945 0.0733 

A) 1.0000 -0.0921 0.1380 

/3i 1.0000 -0.1706 

P2 1.0000 



Table 4. Baycsian parameter estimates for solar cycle 19 using the mean cycle 
as a prior and including data for the entire cycle. 

a be ft Pq /§i 02 

[day" 3 ] [day] [day" 1 ] [day" 1 ] [day- 1 ] [day" 1 ] 

13.23xlO" 8 1401 0.6929 0.0766 18.74 1.891 3.751xl0~ 7 
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Table 5. Bayesian parameter estimates for solar cycle 20 using the mean cycle as 
a prior and including data for the entire cycle. 

a b c k fio j3i /3 2 

[day- 3 ] [day] [day" 1 ] [day" 1 ] [day" 1 ] [day- 1 ] 

8.2605xl0~ s 1400.5 0.8894 0.05624 40.219 1.0995 1.621xl0~ 3 



Table 6. Bayesian parameter estimates for solar cycle 24, using the mean cycle 
as a prior and including all sunspot data from 1 January 2009 to 31 March 2011. 

a b c k fio /9i $2 

[day- 3 ] [day] [day" 1 ] [day- 1 ] [day- 1 ] [day- 1 ] 



4.2962xl0~ 8 1400.0 0.7804 0.09514 10.487 1.6496 0.0040 
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